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Abstract 

We study the long term evolution of the distance between two Kep- 
lerian confocal trajectories in the framework of the averaged restricted 
3-body problem. The bodies may represent the Sun, a solar system 
planet and an asteroid. The secular evolution of the orbital elements 
of the asteroid is computed by averaging the equations of motion over 
the mean anomalies of the asteroid and the planet. When an orbit 
crossing with the planet occurs the averaged equations become singu- 
lar. However, it is possible to define piecewise differentiable solutions 
by extending the averaged vector field beyond the singularity from 
both sides of the orbit crossing set [8], [5]. In this paper we improve 
the previous results, concerning in particular the singularity extrac- 
tion technique, and show that the extended vector fields arc Lipschitz- 
continuous. Moreover, we consider the distance between the Kcplerian 
trajectories of the small body and of the planet. Apart from excep- 
tional cases, we can select a sign for this distance so that it becomes 
an analytic map of the orbital elements near to crossing configurations 
[11]. We prove that the evolution of the 'signed' distance along the 
averaged vector field is more regular than that of the elements in a 
neighborhood of crossing times. A comparison between averaged and 
non-averaged evolutions and an application of these results are shown 
using orbits of near-Earth asteroids. 



1 Introduction 

The distance between the trajectories of an asteroid (orbiting around the 
Sun) and our planet gives a first indication in the search for possible Earth 
impactors. We call it orbit distance and denote it by d m i n } A necessary 

It is often called MOID (Minimum Orbit Intersection Distance, [3]) by the astronomers. 
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condition to have a very close approach or an impact with the Earth is 
that d m i n is small. Provided close approaches with the planets are avoided, 
the perturbations caused by the Earth make the asteroid trajectory change 
slowly with time. Moreover, the perturbations of the other planets produce 
small changes in both trajectories. The value of the semimajor axis of both 
is kept constant up to the first order in the small parameters (the ratio of the 
mass of each perturbing planet to the mass of the Sun). All these effects are 
responsible of a variation of d m i n . We can study the evolution of the asteroid 
in the framework of the restricted 3-body problem: Sun, planet, asteroid. 
Then it is easy to include more than one perturbing planet in the model, in 
fact the potential energy can be written as sum of terms each depending on 
one planet only. 

If the asteroid has a close encounter with some planet, the perturbation 
of the latter generically produces a change in the semimajor axis of the 
asteroid. This can be estimated, and depends on the mass of the planet, 
the unperturbed planetocentric velocity of the small body and the impact 
parameter, see [18]. 

The orbits of near-Earth asteroids (NEAs, i.e. with perihelion distance 
< 1.3 au) 2 are chaotic, with short Lyapounov times (see [19]), at most a 
few decades. After that period has elapsed, an orbit computed by numerical 
integration and the true orbit of the asteroid are practically unrelated and 
we can not make reliable predictions on the position of the asteroid. For this 
reason the averaging principle is applied to the equations of motion: it gives 
the average of the possible evolutions, which is useful in a statistical sense. 
However, the dynamical evolution often forces the trajectory of a NEA to 
cross that of the Earth. This produces a singularity in the averaged equa- 
tions, where we take into account every possible position on the trajectories, 
including the collision configurations. 

The problem of averaging on planet crossing orbits has been studied 
in [8] for planets on circular coplanar orbits and then generalized in [5] 
including nonzero eccentricities and inclinations of the planets. The work in 
[8] has been used to define proper elements for NEAs, that are integrals of an 
approximated problem, see [9]. In this paper we compute the main singular 
term by developing the distance between two points, one on the orbit of 
the Earth and the other on that of the asteroid, at its minimum points. 
This choice improves the results in [8], [5], where a development at the 
mutual nodes was used, because it avoids the artificial singularity occurring 
for vanishing mutual inclination of the two orbits. Moreover, we show that 
the averaged vector field admits two Lipschitz-continuous extensions from 
both sides of the orbit crossing set (see Theorem 4.2), which is useful for the 
numerical computation of the solutions. 

The orbit distance d m % n is a singular function of the (osculating) orbital 

2 1 au (astronomical unit) « 149,597,870 Km 
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elements when the trajectories of the Earth and the asteroid intersect. How- 
ever, by suitably choosing a sign for d m i n we obtain a map, denoted by d m i n , 
which is analytic in a neighborhood of most crossing configurations (see [11]). 

Here we prove that, near to crossing configurations, the averaged evo- 
lution of d m i n is more regular than the averaged evolution of the elements, 
which are piecewise differentiable functions of time. 

The paper is organized as follows. Section 2 contains some preliminary 
results on the orbit distance. In Sections 3, 4, 5 we introduce the averaged 
equations, present the results on the singularity extraction method and give 
the definition of the generalized solutions, which go beyond crossing singu- 
larities. In Section 6 we prove the regularity of the secular evolution of the 
orbit distance. Section 7 is devoted to numerical experiments: we describe 
the algorithm for the computation of the generalized solutions and compare 
the averaged evolution with the solutions of the full equations of motion. We 
also show how this theory can be applied to estimate Earth crossing times 
for NEAs. 

2 The orbit distance 

Let (Ej,Vj), j = 1,2 be two sets of orbital elements of two celestial bodies 
on confocal Keplerian orbits. Here Ej describes the trajectory of the orbit 
and Vj is a parameter along the trajectory, e.g. the true anomaly. We denote 
by £ = (.Ei, £2) the two-orbit configuration, moreover we set V = (1)1,1)2)- 
In this paper we consider bounded trajectories only. 

Choose a reference frame, with origin in the common focus, and write 
Xj = Xj(Ej, Vj), j = 1,2 for the Cartesian coordinates of the two bodies. 

For a given two-orbit configuration £, we introduce the Keplerian dis- 
tance function d, defined by 

T 2 3^4 d(£, V) = \X\ — X 2 \ , 

where T 2 is the two-dimensional torus and | • | is the Euclidean norm. 

The local minimum points of d can be found by computing all the critical 
points of d 2 . For this purpose in [6], [7], [12], [2] the authors have used 
methods of computational algebra, such as resultants and Grobner's bases, 
which allow us to compute efficiently all the solutions. 

Apart from the case of two concentric coplanar circles, or two overlapping 
ellipses, the function d 2 has finitely many stationary points. There exist 
configurations attaining 4 local minima of d 2 : this is thought to be the 
maximum possible, but a proof is not known yet. A simple computation 
shows that, for non-overlapping trajectories, the number of crossing points 
is at most two, see [7]. 
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Let Vh = Vh(£ ) be a local minimum point of V i- >• <i 2 (£, V). We consider 
the maps 



£ (->■ = d(£, V h ) , £ i-)- d min (£) = m.ind h {£) . 

h 

For each choice of the two-orbit configuration £ , d m i n {£) gives the orbit 
distance. 

The maps dh and d m i n are singular at crossing configurations, and their 
derivatives do not exist. We can deal with this singularity and obtain an- 
alytic maps in a neighborhood of a crossing configuration £ c by properly 
choosing a sign for these maps. We note that dh, d m in can present singulari- 
ties without orbit crossings. The maps dh can have bifurcation singularities, 
so that the number of minimum points of d may change. Therefore the 
maps dh, d m i n are defined only locally. We say that a configuration £ is 
non-degenerate if all the critical points of the Keplerian distance function 
are non-degenerate. If £ is non-degenerate, there exists a neighborhood W 
of £ € M 10 such that the maps dh, restricted to W, do not have bifurcations. 
On the other hand, the map d m i n can lose regularity when two local minima 
exchange their role as absolute minimum. There are no additional singular- 
ities apart from those mentioned above. The behavior of the maps dh, d m i n 
has been investigated in [11]. However, a detailed analysis of the occurrence 
of bifurcations of stationary points and exchange of minima is still lacking. 

Here we summarize the procedure to deal with the crossing singularity 
of dh', the procedure for d m % n is the same. We consider the points on the two 
orbits corresponding to the local minimum points Vh = (^i , v^" ) of d 2 : 

x[ h) =X 1 (E 1 , v[ h) ) ; X 2 {h) =X 2 (E 2 , vf ] ) . 

We introduce the vectors tangent to the trajectories E±,E 2 at these points 
(h)_dXi , h) Jh)_dX 2 



and their cross product 
We also define 



-CO _ T C0 x J h ) 

3 —T x X T 2 



A = X x - X 2 , A h = X[ h) - X 2 {h) . 

The vector joins the points attaining a local minimum of d 2 and |Aft| = 
dh- 

From the definition of critical points of d 2 both the vectors t{ , are 
orthogonal to Aft, so that 73 and Aft are parallel, see Figure 1. Denoting 
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Figure 1: Geometric properties of the critical points of d 2 and regularization 
rule. 



by fg , Ah the corresponding unit vectors and by a dot the Euclidean scalar 
product, the distance with sign 

d h = {4 h) ■ A h ) d h (1) 

is an analytic function in a neighborhood of most crossing configurations. 
Indeed, this smoothing procedure fails at crossing configurations such that 
t[ , are parallel. A detailed proof can be found in [11]. Note that, 
to obtain regularity in a neighborhood of a crossing configuration, we lose 
continuity at the configurations with r{ x = and dh 7^ 0. 

The derivatives of dh with respect to each component k = 1 . . . 10 of 
£ are given by 

We shall call (signed) orbit distance the map <2 m j n . 



3 Averaged equations 

Let us consider a restricted 3-body problem with the Sun, the Earth and an 
asteroid. The motion of the 2-body system Sun-Earth is a known function 
of time. We denote by X, X' £ M 3 the heliocentric position of the asteroid 
and the planet respectively. The equations of motion for the asteroid are 



2 2 — 



= — k — -g + fik 



\X'-X\ 3 \X\ 3 



(3) 



where fc is Gauss' constant and fj, is a small parameter representing the ratio 
of the Earth mass to the mass of the Sun. 

We study the motion using Delaunay's elements y = (L,G, Z,£, g, z), 
defined by 

L = k^fa, £ = n(t-t ), 

G = k ^a{\ - e 2 ) , g = u, 

Z = k^J a(l — e 2 ) cos I , z = ft , 
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where (a, e, I, u, £l,£) are Keplerian elements, n is the mean motion and to 
is the time of passage at perihelion. Delaunay's elements of the Earth are 
denoted by (£,', G', Z', g , g', z'). We write £ = (E, E') for the two-orbit 
configuration, where E, E' are Delaunay's elements of the asteroid and the 
Earth respectively. Using the canonical variables y, equations (3) can be 
written in Hamiltonian form as 



where we use 



y = hVyH, 

CO —I 
I o 



(4) 



n G N, 



for the symplectic identity matrix of order In. The Hamiltonian 



H = Hn 



R 



is the difference of the two-body (asteroid, Sun) Hamiltonian 

k 4 



and the perturbing function 



R = nfc 



2L 2 



1 XX' 



X-X'\ |A"I 3 



with X,X' considered as functions of y,y' . 

The function R is the sum of two terms: the first is the direct part of 
the perturbation, due to the attraction of the Earth and singular at colli- 
sions with it. The second is called indirect perturbation, and is due to the 
attraction of the Sun on the Earth. 

We can reduce the number of degrees of freedom of (4) by averaging over 
the fast angular variables £,£', which are the mean anomalies of the asteroid 
and the Earth. As a consequence, t becomes a cyclic variable, so that the 
semimajor axis o is constant in this simplified dynamics. For a full account 
on averaging methods in Celestial Mechanics see [1]. 

The averaged equations of motion for the asteroid are given by 



Y = -h VyR , 



(5) 



where Y = (G,Z,g,zY, Y = {G,Z,g,z) 1 are some of Delaunay's elements, 
and 

W#= TTT^Tn / VyRdldl', 
(2-k) 2 J T 2 
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with T 2 = {(£,£') : —it < £ < n, —it < £' < it}, is the vector of the averaged 
partial derivatives of the perturbing function R. Equation (5) corresponds 
to the scalar equations 

-h_dR -h_dR ^__dR ^__dR 

We can easily include more planets in the model. In this case the per- 
turbing function is sum of terms Ri, each depending on the coordinates of 
the asteroid and the planet i only, with a small parameter /ij, representing 
the ratio of the mass of planet i to the mass of the Sun. 

Note that, if there are mean motion resonances of low order with the 
planets, then the solutions of the averaged equations (5) may be not repre- 
sentative of the behavior of the corresponding components in the solutions 
of (4). 

Moreover, when the planets are assumed to move on circular coplanar 
orbits we obtain an integrable problem. In fact the semimajor axis a, the 
component Z of the angular momentum orthogonal to the invariable plane 3 
and the averaged Hamiltonian H are first integrals generically independent 
and in involution (i.e. with vanishing Poisson's brackets). Taking into ac- 
count the eccentricity and the inclination of the planets the problem is not 
integrable any more. 

In [14] the secular evolution of high eccentricity and inclination asteroids 
is studied in a restricted 3-body problem, with Jupiter on a circular or- 
bit. Nevertheless, crossings with the perturbing planet are excluded in that 
work. In [15] there is a similar secular theory for a satellite of the Earth. 
The dynamical behavior described in [14], [15] is usually called Lidov-Kozai 
mechanism in the literature and an explicit solution to the related equations 
is given in [13]. 

If no orbit crossing occurs, by the theorem of differentiation under the 
integral sign the averaged equations of motion (5) are equal to Hamilton's 
equations 

F = -h VyR (6) 



where 



R= JL / Rdidi' ^ f 



(2tt) 2 Jj2 (2tt) 2 Jj2 \X - X'\ 

is the averaged perturbing function. The average of the indirect term of R 
is zero. 

When the orbit of the asteroid crosses that of the Earth a singularity 
appears in (5), corresponding to the existence of a collision for particular 
values of the mean anomalies. We study this singularity to define generalized 
solutions of (5) which go beyond planet crossings. Since the semimajor axis 



3 Here we mean the common plane of the planetary trajectories. 
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of the asteroid is constant in the averaged dynamics, we expect that the 
generalized solutions can be reliable only if there are no close approaches 
with the planet in the dynamics of equations (4). 



4 Extraction of the singularity 

In the following we denote by £ c a non-degenerate crossing configuration 
with only one crossing point, and we choose the minimum point index h 
such that dh(£ c ) = 0. For each £ in a neighborhood of £ c we consider 
Taylor's development of V \-t d 2 (£,V), V = (£,£')*, in a neighborhood of 
the local minimum point Vh = Vh{£)' 



d 2 (£, V) = d 2 h {£) + -(V - V h ) ■ H h (£)(V - V h ) + Kf\£, V) , 



(7) 



where 



n h (£) 



d 2 d 2 
W 2 



(£,v h (£)) 



is the Hessian matrix of d 2 in Vh = (ih,^) 1 , and 

Tlf\£,V) = Y,ri h \£,V)(V-V h ) a , 



\ a \=a 
3 



^ [\l-t) 2 D a d 2 (£,V h + t(V-V h ))dt 
Jo 



(8) 



(9) 



is Taylor's remainder in the integral form. We introduce the approximated 
distance 

S h = 



< + {y- V h ) ■ A h (V - V h ) 



(10) 



where 



A h = -n h 



N 2 + 



d 2 x 

dt 2 



K\ 2 -^r(E'J' h )-A h 



and 



r h = — (E,£ h ), T> h = —(E',£> h ). 



4 In (8), (9) a = (ai,a 2 ) € (NU {0}) 2 is a multi-index, hence 

|a|=ai+Q2, a\ — ai!«2! , V° = d" 1 ^ 2 , D a f — 
for a vector V = (vi, V2) and a smooth function V i-¥ f(V). 



a |a| / 
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Remark 1. If the matrix Ah is non-degenerate, then it is positive definite 
because Vh is a minimum point, and this property holds in a suitably chosen 
neighborhood W of £ c . At a crossing configuration £ = £ c the matrix «4/j is 
degenerate if and only if the tangent vectors r^, r' h are parallel (see [11]): 

det A h (S c ) = 7fc(£ c ) x r£(£ c ) = . 

First we estimate the remainder function 1/d — \/&h- To this aim we 
need the following: 

Lemma 4.1. There exist positive constants C\, C2 and a neighborhood li 
of (£ c , 

Vfi(£ c )) such that 

C^l <d 2 < C 2 8 2 h (11) 

holds for (£,V) inU. Moreover, there exist positive constants C3, C4 and a 
neighborhood W of £ c such that 

d\ + C 3 \V - V h \ 2 <S 2 h <d 2 h + C A \V - V h \ 2 (12) 

holds for £ inW and for every V £ T 2 . 

Proof. From (8), (9) we obtain the existence of a neighborhood U of (£ c , Vh{£ c )) 
and a constant C5 > such that 

\TZi h) (£,V)\< ]T \ri h \£,V)\\V-V h \ a <C 5 \V-V h \ 3 . (13) 

\a\=3 

We select U so that no bifurcations of stationary points of d 2 occur and 
there exists a constant Cq > with e4(£) > Cq, k / h for each (£, V) G ZY. 
Relation (13) together with (7), (10) yield (11) for some Ci,C 2 > 0. 

Moreover, we can find a neighborhood W of £ c such that there are no 
bifurcations of stationary points of d 2 , and the inequalities (12) hold for 
some Cs^Ci > 0: in fact Ah depends continuously on £ and Ah(£ c ) is 
positive definite. □ 

Proposition 1. There exist C > and a neighborhood W of £ c such that 

<C V (£, V) £ (W x T 2 ) \ Uy, , 

where Uj: = {{£, V h {£)) : £ G £} mf/i S = {«S G W : 4(5) = 0}. 

Proof. By Lemma 4.1 we can choose two neighborhoods W, V of <S C and 
V h (£ c ) respectively such that both (11) and (12) hold inW = WxV. We 



1 1 

d 6 h 
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Figure 2: Sketch for the selection of the neighborhood U = W x V. Here 
Tj = {(£, Vj(£)) : dj{£) = 0} for j = h,k. In this case we restrict W to a 
smaller set (the inner circle), as explained in the proof of Proposition 1. 



restrict W, if necessary, so that there exists Cj > with d > Cj for each 
(£, V) G W x (T 2 \ V) (see Figure 2). \nU\Uj: we have 



1 1 

d S h 



\*l-d 2 \ 



< 



S h d[6 h + d] ~ VC~i[l + VCl] ^ 



for a constant C > 0. Using the boundedness of 1/d, l/5h hi W x (T 2 \ V) 
we conclude the proof. □ 

Now we estimate the derivatives of the remainder function 1/d — 1/dh ■ 

Proposition 2. There exist C > and a neighborhood W of £ c such that, 
if Uk is a component of Delaunay's elements Y, the estimate 



d f\ l_ 

dy k \d S h 



< 



C 



d h +V-V h 



(14) 



holds for each (£,V) £ (WxT 2 )\W E , withU^ as in Proposition 1. Therefore 
the map 

[ JL(l-±-)d£d£', (15) 
J T 2 dy k \d d h J 

where S = {£ € W : dh{£) = 0}, can be extended continuously to the whole 
set W. 



Proof. In the following we denote by Cj,j = 8 . . . 14 some positive constants. 
We write 



9 ( 1 l \ - l ( 1 1 \ d5 h 1 dn 



(/') 



<9y fc Vd <^V 2 d 3 Jdy k 2d 3 dy A 
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and give an estimate for the two terms at the right hand side. We choose a 
neighborhood U = W x V of (£ c , Vh{£ c )) as in Proposition 1 so that, using 
(11), (12) and the boundedness of the remainder function, we have 



1 


1 






1 












U \ Uy, ■ Moreover 




dsl 


< 


ddl 




dyk 




dyk 



1 



1 + i)< 



c 8 



6 2 h + S h d ' <P) " d 2 h + \V-V h \ 2 



+ C 9 \V-V h \<C w (d h + \V-V h \) 



since 



2^ ■ A h (V - V h ) + (V- V h ) ■ ^(V - V h ) , (16) 
dyk dy k dy k dy k 



and the derivatives 
dV k 



are uniformly bounded for £ £ W since bifurcations do not occur. 
Hence the relation 











d 3 ; 





2Cu 



4 + l^-^l 2 " 4 + l^-Vfc| 



(17) 



holds in with Cn = CsCio- We also have 



an 



(h) 



dyk 



<C l3 \V-V h \ 2 , 



(18) 



for (£, V) G Uq, in fact 



sup \r£'\ < +oo , sup 
u u 



dyk 



< +oo , 



(19) 



for each a = (01,02) with \a\ = 3. Using again (11), (12) we obtain 



1 dll 



00 



d 3 dy k 



< 



d h + \V-V h \ 



(20) 



From (17), (20) we obtain (14) and the assert of the proposition follows using 
the boundedness of ^(l/d) , ^(l/5 h ) in W x (T 2 \ V). □ 
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Prom (14) in Proposition 2 the average over T 2 of the derivatives of 1/d— 
l/5h in (15) is finite for each £ in W, and can be computed by exchanging the 
integral and differential operators: therefore the average of the remainder 
function is continuously differentiable in W. 

On the other hand, the average over T 2 of the derivatives with respect 
to yk of l/b~h are non-convergent integrals for £ € S: for this reason the 
averaged vector field in (5) is not defined at orbit crossings. Next we show, 
exchanging again the integral and differential operators, that the average of 
these derivatives admit two analytic extensions to the whole W from both 
sides of the singular set S. 

For this purpose, given a neighborhood W of £ c , we set 

w+ = w n {d h > 0} , r = wn {4 < 0} , 

with dh given by (1). 

Proposition 3. There exists a neighborhood W of £ c such that the maps 

W + 3£^^-[ ^-d£d£', W~3£^^- f ^d£d£', 
oyk Jt 2 Oh oy k J T 2 d h 

where y k is a component of Delaunay's elements Y, can be extended to two 
different analytic maps Qiu,QZu, defined in W. 

Proof. We choose W as in Proposition 2 and, if necessary, we restrict this 
neighborhood by requiring that r} x 7^ in W, so that d~h\w is analytic. 
To investigate the behavior close to the singularity, for each £ S W, we can 
restrict the integrals to the domain 

V = V(V h ,r) = {V e T 2 : (V - V h ) ■ A h (V - V h ) < r 2 } , (21) 

1/2 

for a suitable r > 0. By using the coordinate change £ = A h (V — Vh) and 
then polar coordinates (p,0), defined by (pcos 6, psin#) = £, we have 

/ -j- dldJt = j=t A [ . 1 di 
Jv Oh ^det A h J b / d 2 + |^|2 

with £> = {££ R 2 : |£| < r}. The term -2ird h /y/dei A h is not differentiable 
at £ = £ c £ S. We set 

d / 2tt \ / 2 2 27T <4 g4 

^' fc "^AVdeT^y , V^ +r + v^it^y^T^^ 
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with dh as in (1), and define on W the two analytic maps 



'h,k 







(- 



2vr 



dy k V Vdet At 



4=F 



2vr 3d 



h 







Vdet A h dy k dy k J T 2\ V S h 



1 



+ t— / ^dldl' . (22) 



We observe that Q^ k (resp. Q h k ) corresponds to the derivative of fj 2 l/Shdl di 
with respect to y k on W + (resp. W~). □ 

Now we state the main result. 

Theorem 4.2. The averages overT 2 of the derivatives of R with respect to 
Delaunay's elements y k can be extended to two Lipschitz-continuous maps 
{§^) h on a neighborhood W of £ c . These maps, restricted to W + , W~ 

respectively, correspond to . Moreover the following relations hold: 



Diff, 



dR \ de_f (dRy (dR\+ _ 
dyk) \dy k )h \dy k )h 

yk 2 \ d ( 1 
7T [dy k \^/det(A h ) 

with the derivatives of dh given by (2). 



d h + 



dd h 



^det(A) d Vk 



, (23) 



Proof. Using the results of Propositions 2, 3 we define the extended maps 
by 



( 



dR\± 



fik 2 



\dy k Jh (2vr) 5 



d 

'2 dy k 



1 1 

d S h 



* + Qik 



with Qh k given by (22). We show that the maps £ i— >■ (J^)^ 1 ?) are 

Lipschitz-continuous extensions to W of J^-. The maps Q^ k are analytic in 

W, thus we only have to study the integrals f j2 df^(V^ — V^O d£d£'. From 
Proposition 2 we know that these maps are continuous. 

We only need to investigate the behavior close to the singularity, therefore 
we restrict these integrals to the domain V introduced in (21). We prove 
that the maps 



d 1 

dy k S h 



didl 1 , 



/ -^--dtdJt!, 
oyj Jv oy k d 



W\^3£^^- [ 
dyj Jx 

with j = 1 ... 4, are bounded. First observe that the derivatives 
d 1 



_d_ 

dyj Jv dyk Sh 



did? 



3 1 851 < 



11 dHl 



v\^Sl dyj dy k 2 S\ dyjdy k 



did? 



(24) 
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are bounded in W \ S, otherwise we could not find the analytic extensions 
QtuiQ^k introduced in Proposition 3. 5 Then we show that the maps 







d 1 



dyj Jv dyk d 



dtdl' 



3 1 dd 2 dd 2 1 1 d 2 d 2 



v \4:d 5 dyj dy k 2 d 3 dyj dy k 



did? (25) 



are bounded in W \ S. Using (7), (10) we write the integrand function in 
the right hand side of (25) as the sum of 



4 5 5 h dyj dy k 
and of terms of the following kind: 



1 i d 2 5 2 

2 Si dyjdy k 



(26) 



3 1 an 



(h) 



4 Si d Vj 



dU 



(h) 



dyk 



1 + P W) +P W°& 
5 ' 5 dy k 



1 V { 3 h) d 2 K 3 

2 Si dyjdy k 



(fc) 



(27) 



3 1 dK 



(h) 



A5\ dyj dy k 
4 Si dyj dy k ' 



1 1 d 2 iz[ h) 

2 5\ dyjdy k 

17f } d 2 S 2 
2 5\ dyjdy k 



(28) 



(29) 



The integrals over T> of the terms in (26) are not bounded in W \ E, but 
their sum is bounded and corresponds to (24). In the following we denote 
by Cj,j = 15 . . . 34 some positive constants. Moreover we use the relation 
d 2 = 5 2 + 1Z 3 and the developments 



(5i+n { 3 h) y/ 2 $ s h l 



i + vi h) 



(5 = 3,5) 



with 



v w = r ( h){£j v)= ^2 P f>(s, v)(v - v h y 

101=1 

00, 

(f,^i + t(F-Vft))d< 



1 + — ^ 



si 



3 Actually we can prove that 



(30) 



where H^,^' 

are bounded in W \ E, and 



rrO) _ 



'2tt 



dd h dd h dV h dV h 



is unbounded but cancels out in the difference. 
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By developing (30) we obtain 



1 + 



K 



(/<) 



SI 



n 



(h) 



V si 



{£,V h + t(V-V h ))dt 



13 '^3 



{£, V h + t(V - V h )) dt + SR^(£ , V) 



with V)| < C 1S \V -V h \,s = 3, 5. Moreover, we have 



SI 



We can estimate the terms in (31) as follows: 



CO 



l«|=3 



where 
so that 
Moreover 
so that 

We conclude that 
D@ I 



D fi r^\ < C 16 , |l/(U - V h ) a \ < C 17 \V - V h \ 2 , 



\D p nf ] \ <c 18 \v-v h \ 2 



dH\ = 1D\V - V h ) ■ A h {V - V h ) 

\dH h \ <c 19 \v-v h \ . 



K. 



{hy 



si 



<C 20 , so that \pf^£,V)\<C 21 , 



and we obtain the estimate 

\vi h \£,V)\<C 22 \V-V h \ 
for (£,V) £ U Q . Using (16), (18), (32) and the estimate 



d 2 n{ h) 



dyjdy k 



<c 23 \v-v h \, 



that follows from the boundedness of 

r (h) 



dri h) d 2 A h) 



dv h d 2 v h 



dyk ' dyjdyk' dy k ' dyjdy k 
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(31) 



(32) 



we can bound both terms in (27) by C24/IV — Vh\, which has finite integral 
over P. 6 

To estimate the integrals of the terms in (28) we observe that 



dK 



th) 



9 Vj 
with 



/ . r a,0 
|a|=3 



dyjdy k 



E> 

|o|=3 



Sh) d 2 (v-v h y 



a,0 



dyjdy k 



+6 



(h) 
2 > 



r£J = ) = rW(£, ^(f )) , | 6 < h) | < C 25 |V - K fc |« (• = 2, 3) . 



Then, using (16) and writing dV for dldt' , we have 



1 dKi h) 051 dv 



v 5\ dyj dy k 



< 



ddl 



dy k 



V \r (h) \ 

/ ; \ r a,0\ 
\a\=3 



V 



1 d(v - v h y 



dV + 



T> S h 



+2 £ K 

|a|=3 



00 I 

ol 



1 a(y - v h ) 



2? °/i 



A(^ - V h 



dV 


+ c 2& [ 




Jv 



\V-Vh\ l 



(33) 



dV 



and 



1 d 2 ^ 

■d 8\ dyjdy k 



dV 



< E kJSl 

|a|=3 



1 d 2 (^ - Vfe 



+C27 



IV — T4 12 



<5? 



■dV. 



1/2 

Passing to polar coordinates (p,0), defined by (p cos 9, p sin 0) = .4^' (V 
Vft), we find that 



(34) 



94 




1 d(v - v h ) a 




X 





for each £ G W \ E and a with |a| =3, in fact 



54 




dyk 


L 



16 



W 1 



K + p2 )5 /2 



dp<^° (* = 3,4) 
«/). 



(35) 



Moreover, passing to polar coordinates (p,0), we have 



p 



1 d(v - v h y 

tfjj d yj 

A 



dVh 

dyk 



A h (V - V h 



dV 



/ (A2 I 2^5/2 ^ COS 6 ^ ^ Sin ^ ^ = ^ 

JO \ a h + P ) | 7 | =3 -/0 



6 The boundedness of ^ ^ on W follows by differentiating with respect to yj the 
relation 



oy k dy k 
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for some functions c 7 : W \ £ — > K, 7 = (71,72)- Thus the integrals in (33) 
are uniformly bounded in W \ X. In (36) we have used 



2- 



(cos <9) 71 (sin 6y /2 d6 = 



(37) 



for each 7, with odd I7I =71+72 • Finally, using again (37), we obtain 



1 d 2 (v-v h y 

v Sl dyjdy k 



dV 



< 



E & , 

l7l=l 



2ti- 



{cos ey 1 (sm ey 2 d9 



p 



+ Q 



31 



1 



V 



\v-v h \ 



dv = a 



31 



1 



V 



\v - v h \ 



(dl + p 2 f/ 2 

dV 



dp 



for some functions 6 7 : W \ £ — > K. Hence also the integrals in (34) are 
uniformly bounded in W \ X. 

To estimate the integrals of the terms in (29) we make the following decom- 
position: 

Jh) _ (h) , (h) 



Iks + 



,s ' 



with 



a (h) 

s 

' ~ 2 



£ 



,.00 
a,0 



\a\=3 



1 / D> 3 1) 2 

- D P(v-v h r-—^(v-v h y 



si 



(S,V h + t{V -V h ))dt 



and \w$\ < C32I — Vh\- For the first term in (29) we obtain 



ri h) 351 K dy 



f K 

Jv SI 



ddl dd\ 



+ 4 



dyj dy k 



101=1 



dyj dy k 
dVh 



jv h \p\=i 



A h (V - V h ) 



dV h 
dyk 



A h {V - V h ) 



(38) 
dV 



+ C33 



where we have used polar coordinates and the inequalities (35). The two 
integrals at the right hand side of (38) vanish: in fact using Fubini-Tonelli's 
theorem and passing to polar coordinates (p,6), by relations (37) we obtain 

= E E 

|/3|=l| 7 |e{3,5}' 



1 t-r 



JO 



2tt 



(p, t) dp dt / (cos #) 71 (sin #) 72 d6 = 



17 



for some functions <^g~ : M + x R — >• R. The computation for the other 
integral is analogous. 

The second term in (29) is estimated in a similar way: 




and the integral at the right hand side vanishes as well. 

We conclude the proof observing that, using (22) and the theorem of 
differentiation under the integral sign, the derivatives {§^ L ) h : {.W^)h^ re ~ 

stricted to W + , respectively, correspond to J^-, and their difference in 
W is given by (23). □ 

Remark 2. If £ c is an orbit configuration with two crossings, assuming that 
dh{£c) = for h = 1,2, we can extract the singularity by considering the 
approximated distances <5i , 82 and the remainder function \/d—\/8\ — \ /S2 ■ 



5 Generalized solutions 

We show that generically we can uniquely extend the solutions of (5) beyond 
the crossing singularity d m i n = 0. This is obtained by patching together 
classical solutions defined in the domains W + with solutions defined in W~ , 
or vice versa. 

Let a > be a value for the semimajor axis of the asteroid and Y : I — > M 4 
be a continuous function defined in an open interval Jcl, representing a 
possible evolution of the asteroid orbital elements Y = (G,Z,g,z). We 
introduce 

£(t) = (E(t),E'(t)), (39) 

with 

E(t) = (kVZ,Y(t)) , (40) 

where k is Gauss' constant and e' is a known function of time representing 
the evolution of the Earth. 7 

Let T(Y) be the set of times t c € I such that d m i n (£(t c )) = 0, and 
assume that each t c is isolated, so that we can represent the set 

/ \ T(Y) = Ujsjflj 

7 In the case of one perturbing planet E' (t) is constant and represents the trajectory of 
a solution of the 2-body problem. If we consider more than one perturbing planet then 
E (t) changes with time due to the planetary perturbations. 
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as disjoint union of open intervals Ij, with M a countable (possibly finite) 
set. 

Definition 5.1. We say that Y is a generalized solution of (5) if its restric- 
tion to each Ij,j G TV is a classical solution of (5) and, for each t c G T(Y), 
there exist finite values of 

lim Y(t) , lim Y{t) . 

t-¥tt t-¥tc 

Choose Yq G M 4 and a time to such that d m i n (£o) > 0, with £q = (Eq, E' ), 
Eq = (ky/a,Yo), E'q = E'(to). We show how we can construct a generalized 
solution of the Cauchy problem 

Y=-hWR, Y(t Q ) = Y . (41) 

Let Y(t) be the maximal classical solution of (41), defined in the maximal 
interval J. Assume that t c = sup J < +oo, and lim. .- £ (t) = £ c , with 
£ c a non-degenerate crossing configuration such that d m i n {£ c ) = dh{£ c ) = 
for some h. Let W, W 1 * 1 be chosen as in Theorem 4.2. Suppose that there 
exists r G (t ,t c ) such that £(t) G W + for t G (r,t c ). Let Y T = F(r). By 
Theorem 4.2 there exists Y c G M such that 

lim Y(t) = Y c . (42) 

t-ttc 

In fact relation (42) is fulfilled by the solution of the Cauchy problem 8 



Y = -J 2 (Vyi?)+ , Y{t) = Y T , (43) 

which corresponds to the solution of (41) in the interval (r,t c ) and is defined 
also at the crossing time t c . Let us denote by Y c its value for t = t c . Using 
again Theorem 4.2 we can extend Y(t) beyond the crossing singularity by 
considering the new problem 

t = -J 2 (WR)h , Y(t e ) = Y c . (44) 
The solution of (44) fulfils 

lim F(t) = Y c + Dift h (WR)(£ (Q) ■ (45) 
t->t+ 

The vector field in (44) corresponds to — J2Vyi? on W~, thus we can continue 
the solution outside W and this procedure can be repeated at almost every 

8 Here (VyR)^ is the vector with components (J^)^ , k — 1 . . .4 introduced in Theo- 
rem 4.2. 
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crossing singularities. Indeed, the generalized solution is unique provided 
the evolution t \— >■ £ (t) is not tangent to the orbit crossing set X. 
Moreover, if det^A^ = the extraction of the singularity, described in Sec- 
tion 4, cannot be performed. 

In case £(t) G W~ for t G (r, t c ) the previous discussion still holds if we 
exchange (Vyi?)^ with (Vy-R)7. In this case (45) becomes 

lim t(t) = Y C - Dm h (WR)(S(tc)) ■ 

6 Evolution of the orbit distance 

We prove that the secular evolution of d m i n is more regular than that of the 
orbital elements in a neighborhood of a planet crossing time. We introduce 
the secular evolution of the distances dh and of the orbit distance d m i n : 

d h {t) = d h (E{t)) , d min {t) = d min (£(t)) . (46) 

Assume these maps are defined in an open interval containing a crossing 
time t c , and suppose £ c = £(t c ) is a non-degenerate crossing configuration 
at time t c , as in Section 4. 

In the following we shall discuss only the case of dh- The same result 
holds for d m in, taking care of the possible exchange of role of two local 
minima dh,dk as absolute minimum. 

Proposition 4. Let Y(t) be a generalized solution of (41) and £{t) as in 
(39), (40). Assume t c G T(Y) is a crossing time and £ c = £(t c ) is a non- 
degenerate crossing configuration with only one crossing point. Then there 
exists an interval (t a ,tb), t a < t c < t^ such that dh G C 1 ((t a , tb); K). 

Proof. Let the interval (t a ,tb) be such that £((t a ,tb)) C W , where W is 
chosen as in Theorem 4.2. We can assume that £{t) G W + for t G (t a ,t c ), 
£(t) G W~ for t G (t c ,tb) (the proof for the opposite case is similar). For 
t G (t a ,tb) \ {t c } the time derivative of dh is 

\(t) = V £ d h (£(t)) ■ t(t) = V Y d h (£(t)) • Y(t) + V E ,d h (£(t)) • E(t) 
= - (Vydh ■ hWR) (S(t)) + V E 'd h (£(t)) ■ E(t) . 

Here Vg, Vy, V^y denote the vectors of partial derivatives with respect to 

£, Y, E' respectively. The derivative E (t) exists also for t = t c . On the other 
hand, by Theorem 4.2, the restrictions of Vy /?(£(£)) to t < t c and t > t c 
admit two different continuous extensions to t c . By (23), since dh(£(t c )) = 0, 
we have 

lim dh(t) — lim dh(t) 



Diff h (Vyi?) -hVydh 
fik 2 



£=£ c 



7T\/det Ah 



{d h ,dh}y 



e=£ c 
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where {, }y is the Poisson bracket with respect to Y. Thus the time deriva- 
tive of dh exists and is continuous also in t = t c . □ 
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7 Numerical experiments 



7.1 The secular evolution program 

Using a model with 5 planets, from Venus to Saturn, we compute a planetary 
ephemerides database for a time span of 50, 000 yrs starting from epoch 
MJD (November 17, 1858) with a time step of 20 yrs. The computation is 
performed using the FORTRAN program orbit9, included in the OrbFit 
free software 9 . From this database we can obtain, by linear interpolation, 
the evolution of the planetary trajectories at any time in the specified time 
interval. 

We describe the algorithm to compute the solutions of the averaged equa- 
tions (5) beyond the singularity, where R is now the sum of the perturbing 
functions R4, i = 1 . . . 5, each related to a different planet. We use a Runge- 
Kutta-Gauss (RKG) method to perform the integration: it evaluates the 
averaged vector field only at intermediate points of the integration time in- 
terval. When the asteroid trajectory is close enough to an orbit crossing, 
then the time step is decreased to reach the crossing condition exactly. 

From Theorem 4.2 we can find two Lipschitz-continuous extensions of 
the averaged vector field from both sides of the singular set S. 

To compute the solution beyond the singularity we use the explicit for- 
mula (23) giving the difference between the two extended vector fields, either 
of which corresponds to the averaged vector field on different sides of E. We 
compute the intermediate values of the extended vector field just after cross- 
ing, then we correct these values by (23) and use them as approximations of 
the averaged vector field in (5) at the intermediate points of the solutions, 
see Figure 3. This RKG algorithm avoids the computation of the extended 
vector field at the singular points, which may be affected by numerical in- 
stability. 

A difficulty in the application of this scheme is to estimate the size of 
a suitable neighborhood W of the crossing configuration £ c fulfilling the 
conditions given in Section 4. 

7.2 Comparison with the solutions of the full equations 

We performed some tests to compare the solutions of the averaged equa- 
tions (5) with the corresponding components of the solutions of the full 
equations (4). Here we show two tests with the asteroids 1979 XB and 
1620 (Geographos). We considered the system composed by an asteroid 
and 5 planets, from Venus to Saturn. We selected the 8 values kn/A, with 
k = ... 7, for the initial mean anomaly of the asteroid and the same for 
the planets. Using the program orbit9, we performed the integration of 
the system with these 64 different initial conditions (i.e. we selected equal 

9 http : / /adams .dm.unipi . it/~orbmaint/ orbf it/ 
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Figure 3: Runge-Kutta-Gauss method and continuation of the solutions of 
equations (5) beyond the singularity. The crosses correspond to the inter- 
mediate values. 
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Figure 4: Averaged and non-averaged evolutions of asteroid 1979 XB. 



initial phases for all the planets). Then we considered the arithmetic mean 
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Figure 5: Averaged and non-averaged evolutions of asteroid 1620 (Ge- 
ographos). 

of the four equinoctial 10 orbital elements h, k,p, q of the asteroid over these 
evolutions, and compared them with the results of the secular evolution. In 
Figures 4, 5, we show the results: the crosses indicate the secular evolution, 
the continuous curve is the mean of full numerical one and the gray region 
represents the standard deviation from the mean. The correspondence be- 
tween the solutions is good. During the evolution the distance between the 
asteroid and the Earth for some initial conditions attains values of the order 
of 10~ 4 au for 1620 (Geographos), and 10~ 3 au for 1979 XB. In Figure 5 the 
Earth crossing singularity is particularly evident near the epoch 3000 AD. 

Some numerical tests of the theory introduced in [8] , with the planets on 
circular coplanar orbits, can be found in [10]. 

10 We recall that 

h — e sin(w + fl) , k = e cos(cj + fi) , p = tan(//2) sin(fi) , q — tan(//2) cos(fi) . 
The equinoctial orbital elements have been introduced in [4]. 
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7.3 An estimate of planet crossing times 

The results of Section 6 can be used to estimate the epoch in which the orbit 
of a near-Earth asteroid will cross that of the Earth. We are interested in 
particular to study the behavior of those asteroids whose orbits will cross 
the Earth in the next few centuries, so that they must have a small value 
of d m i n already at the present epoch. We can consider, for example, the set 
of potentially hazardous asteroids (PHAs), which have d m i n < 0.05 au and 
absolute magnitude H mag < 22, i.e. they are also large. 

In Figure 6 we show 3 different evolutions of the signed orbit distance 
dmin for the PHA 1979 XB. Here we draw the full numerical (solid line), sec- 
ular (dashed) and secular linearized (dotted) evolution of d m i n . By Propo- 
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0.02 

0.015 

S 0.01 
< 

(D 

j§ 0.005 

! o 

-0.005 

-0.01 
-0.015 

-0.02 

1950 2000 2050 2100 2150 2200 

time (yr) 

Figure 6: Different evolutions of d m i n for 1979 XB: full numerical (solid line), 
secular (dashed) and secular linearized (dotted). 

sition 4 the linearization of the secular evolution d m i n (t) can give a good 
approximation also in a neighborhood of a crossing time. 

We propose a method to compute an interval J of possible crossing times. 
We sample the line of variation (LOV), introduced in [16], which is a sort of 
'spine' of the confidence region (see also [17]), and compute the signed orbit 
distance d m i n for each virtual asteroid (VA) of the sample. Then we compute 
the time derivative of d m i n for each VA and extrapolate the crossing times 
by a linear approximation of the evolution. We set J = [ti , ^2] ? with £i,t2 
the minimum and maximum crossing times obtained (see Figure 7). In the 
computation of J we take into account a band centered at the Earth crossing 
line d m i n — 0: in this test the width of the considered band is 2 x 10~ 3 au. 
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Figure 7: Computation of the interval J (horizontal solid line) for asteroid 
1979 XB. The transversal solid line corresponds to the linearized secular 
evolution of the nominal orbit. The linearized secular evolution of the VAs 
are the dotted lines. 

We describe a method to assign a probability of occurrence of crossings 
in a given time interval, which is related to the algorithm described above. 
For each value of the LOV parameter s we have a VA at a time t, so that 
we can compute d m i n (t). Thus, using the scheme of Figure 7 we can define 
a map T from the LOV parameter line to the time line. The map T gives 
the crossing times, using the linearized secular dynamics, for the VAs on the 
LOV that correspond to the selected values of the parameter s. Moreover, 
we have a probability density function p(s) on the LOV. Therefore, given an 
interval I in the time line, we can consider the set Uj = T _1 (I) and define 
the probability of having a crossing in the time interval I as 



Finally, in Figure 8 we show the corresponding interval J' obtained by 
computing the secular evolution (without linearization) of the orbit distance 
for each VA of 1979 XB. The sizes of J and J' are almost equal, but the left 
extremum of J' is ~ 10 years before. 
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Figure 8: Computation of the interval J' (horizontal solid line) for asteroid 
1979 XB. The enhanced transversal curves refer to the nominal orbit: solid 
line corresponds to secular evolution, linearized is dashed. The dotted curves 
represent the secular evolution of the VAs. 

8 Conclusions and future work 

We have studied the double averaged restricted 3-body problem in case of 
orbit crossing singularities, improving and completing the results in [8], [5]. 
This problem is of interest to study the dynamics of near-Earth asteroids 
from a statistical point of view, going beyond the Lyapounov times of their 
orbits. We have also proved that generically, in a neighborhood of a crossing 
time, the secular evolution of the (signed) orbit distance is more regular than 
the averaged evolution of the orbital elements. 

The solutions of this averaged problem have been computed by a numer- 
ical method and then compared with the solutions of the full equations in 
a few test cases. The results were good enough; however, we expect that 
the averaging technique fails in case of mean motion resonances or close en- 
counters with a planet. We plan to perform numerical experiments with a 
large sample of near-Earth asteroids showing different behaviors: this will 
be useful to understand the applicability of the averaging technique to the 
whole set of NEAs. 



27 



Acknowledgments 

We would like to thank an anonymous referee for his interesting comments, 
that allowed us to improve the paper. 

References 

[1] (MR1656199) V. I. Arnold, V. V. Kozlov and A. I. Neishtadt, "Mathe- 
matical Aspects of Classical and Celestial Mechanics," Springer, 1997. 

[2] (MR2162789) [10.1007/sl0569-004-3207-l] R. V. Baluyev and K. V. 
Kholshevnikov, Distance between Two Arbitrary Unperturbed Orbits, 
Cel. Mech. Dyn. Ast., 91 (2005), 287-300. 

[3] E. Bowell and K. Muinonen, Earth-crossing asteroids and comets: 
Groundbased search strategies, in "Hazards due to Comets & Asteroids" 
(Ed. T. Gehrels), Tucson: The University of Arizona Press, (1994), 149- 
197. 

[4] R. A. Broucke and P. J. Cefola, On the equinoctial orbit elements, Cel. 
Mech. Dyn. Ast., 5 (1972), 303-310. 

[5] (MR1956518) [10.1023/A:1020178613365] G. F. Gronchi, Generalized 
averaging principle and the secular evolution of planet crossing orbits, 
Cel. Mech. Dyn. Ast., 83 (2002), 97-120. 

[6] (MR1924415) [10.1137/S1064827500374170] G. F. Gronchi, On the 
stationary points of the squared distance between two ellipses with a 
common focus, SIAM Journ. Sci. Comp., 24 (2002), 61-80. 

[7] (MR2186811) [10.1007/sl0569-005-1623-5] G. F. Gronchi, An algebraic 
method to compute the critical points of the distance function between 
two Keplerian orbits, Cel. Mech. Dyn. Ast., 93 (2005), 297-332. 

[8] (MR1675692) [10.1023/A:1008315321603] G. F. Gronchi and A. Milani, 
Averaging on Earth-crossing orbits, Cel. Mech. Dyn. Ast., 71 (1998), 
109-136. 

[9] G. F. Gronchi and A. Milani, Proper elements for Earth crossing 
asteroids, Icarus, 152 (2001), 58-69. 

[10] G. F. Gronchi and P. Michel, Secular orbital evolution, proper elements 
and proper frequencies for near-earth asteroids: A comparison between 
semianalytic theory and numerical integrations, Icarus, 152 (2001), 
48-57. 



28 



[11] (MR2291871) [10.3934/dcdsb. 2007.7.755] G. F. Gronchi and G. Tom- 
mei, On the uncertainty of the minimal distance between two confocal 
Keplerian orbits, Discrete Contin. Dyn. Syst. Ser. B, 7 (2007), 755-778. 

[12] (MR1750216) [10. 1023/A: 1008312521428] K. V. Kholshevnikov and N. 
Vassiliev, On the distance function between two keplerian elliptic orbits, 
Cel. Mech. Dyn. Ast., 75 (1999), 75-83. 

[13] (MR2317260) [10.1007/sl0569-007-9069-6] H. Kinoshita and H. Nakai, 
General solution of the Kozai mechanism, Cel. Mech. Dyn. Ast., 98 
(2007), 67-74. 

[14] Y. Kozai, Secular perturbation of asteroids with high inclination and 
eccentricity, Astron. Journ., 67 (1962), 591-598. 

[15] M. L. Lidov, The evolution of orbits of artificial satellites of planets 
under the action of gravitational perturbations of external bodies, Plan. 
Spa. Sci., 9 (1962), 719-759. 

[16] A. Milani, S. R. Chesley, M. E. Sansaturio, G. Tommei and G. Valsecchi, 
Nonlinear impact monitoring: Line of variation searches for impactors, 
Icarus, 173 (2005), 362-384. 

[17] (MR2778686) A. Milani and G. F. Gronchi, "Theory of Orbit Deter- 
mination," Cambridge Univ. Press, 2010. 

[18] G. B. Valsecchi, A. Milani, G. F. Gronchi and S. R. Chesley, Resonant 
returns to close approaches: Analytical theory, Astron. Astrophys., 408 
(2003), 1179-1196. 

[19] A. Whipple, Lyapunov times of the inner asteroids, Icarus, 115 (1995), 
347-353. 



29 



